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Abstract 

Using finite-size scaling techniques, we study the critical properties 
of the site-diluted Ising model in four dimensions . We carry out a 
high statistics Monte Carlo simulation for several values of the dilu- 
tion. The results support the perturbative scenario: there is only the 
Ising fixed point with large logarithmic scaling corrections. We ob- 
tain, using the Perturbative Renormalization Group, functional forms 
for the scaling of several observables that are in agreement with the 
numerical data. 
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1 Introduction 



A possible way to obtain new Universality Classes (UC) is to add disorder to 
known pure systems. The Harris criterion says that, if the specific heat 
diverges with a power law in the pure system, then, the disorder will change 
the critical behavior of the model, i.e. a new UC will appear. Conversely, 
if the specific heat does not diverge in the pure system, then, the critical 
exponents of the disordered system will remain unchanged. In the limiting 
case, what amounts to a discontinuous or a logarithmically divergent specific 
heat, the criterion does not apply and we need to study analytically or/and 
numerically the system. 

The quest and the characterization of new UC are very important in 
dimensions two and three (with a direct relevance on Condensed Matter 
Physics) and in four dimensions (with implications on High Energy Physics) . 
In the last case it is crucial to characterize all the possible UC, in order to be 
able to define a Field Theory on a non-perturbative basis. As the Gaussian 
model gives a trivial one (i.e. at long distances the theory will be free) we 
are interested in finding a non-Gaussian UC. 

In this paper we will study the four dimensional site-diluted Ising model 
that was previously studied [0 by two of the authors, who calculated numeri- 
cally the critical exponents, analyzing the divergences with the temperature. 
Their results pointed to non Gaussian critical exponents, for large values 
of the dilution, but it was noted the possibility of a crossover between the 
behavior found and the Gaussian one. 

In order to obtain accurate measures of the critical properties we have 
repeated the simulations in a greater number of spin configurations. The use 
of Finite Size Scaling (FSS) techniques allows us to work in large lattices at 
the critical point. 

Using the Perturbative Renormalization Group (PRC) equations we cal- 
culate the dependence of the observables at the critical point with the scale 
of the system including logarithmic corrections. These are different from the 
pure system. 

Our determination of the critical exponents and other critical properties 
matches very well with the predictions of the PRC: a Gaussian critical be- 
havior with logarithmic corrections. We check this behavior along the critical 
line in a wide range of concentrations, from p = 0.8 to p = 0.3 (the percola- 
tion threshold is near 0.2). We remark that a scenario based on hyperscaling 
seems completely unlikely from our numerical simulations. 
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The structure of the paper is as follows. In the next section we define the 
model and the observables. In section 3 we present our analytical calcula- 
tions, obtaining the FSS formulas and calculating the values of two different 
Binder cumulants in the thermodynamic limit. In section 4 we describe the 
numerical methods and the different techniques that we will use to analyze 
the observables. In section 5 we show our numerical results confronting them 
with the analytical predictions. Finally we report the conclusions. 

2 The model 

The model we study numerically in this article, is defined in terms of Z2 
spin-variables placed in the nodes of a hypercubic four-dimensional lattice. 
The action is 

S ^ eiejGiaj , (1) 

where the sum is extended over nearest-neighbors, and the are quenched, 
uncorrelated random variables, taking the value 1, with probability p, and 
with probability 1 — p. An actual {e,} configuration, will be called a sample 
from now on. For every observable it is understood that first one performs 
the Ising model calculation and then the e-average. 

In the following, we shall denote an Ising average with brackets, while 
the sample average will be overlined. The observables will be denoted with 
calligraphic letters, i.e. O and with itahcs the double average O — {O). We 
define the total nearest-neighbor energy and the normalized magnetization 
as 

£ = ^ eiGiejaj , M = ^^eiGi , (2) 

V being the volume (defined as L*, where L is the lattice size). We also 
define the susceptibility as 

X = VjA^ ■ (3) 
Another very useful quantity is the Binder parameter: 
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Other kind of Binder parameter, meaningless for the pure system, can be 
defined as 



92 = 2 • 5) 

A very convenient definition of the correlation length in a finite lattice, 
reads 

(-(M^)', (6) 



^4sin^(7r/L) 

where F is defined in terms of the Fourier transform of the magnetization 



Hk) = ^J2^''"''r(rr , (7) 



V 

r 



as 

V 



F = —(|J^^(27r/L,0, 0,0) |2 + permutations) . (8) 

This definition is very well behaved for the finite-size scaling (FSS) method 
we employ 0, and it is also fair natural for considerations about triviality 0. 
Finally, we measure the specific heat 



C = V^\S^) - {S)^ . (9) 



3 Analytical predictions 

The field-theoretical study of the model (|l]) is performed by means of a 0^ 
theory with a random mass term, whose action is 

S[<P] = J d'^x Q {d,<pf + ^m'ix)<P' + 1^;^^^ . (10) 

Here it is assumed that the mass term is a quenched, spatially-uncorrelated, 
stochastic variable. We will later argue that the only relevant parameters of 
the disorder distribution are its mean, r, and variance, zA^, so we will assume 
for simplicity the Gaussian distribution 
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dP[m^(x)] oc d[m^] exp 



(m^(x) 



(11) 



When the disorder is quenched we need to compute, in a first step, the free 
energy of the system for a given choice of the disorder (in this case of the mass 
term), and then average this free energy with the probabihty distribution of 
the disorder. To manage this kind of problems it is very useful to use the 
so-called replica trick p. 

Let us introduce n replicas of the initial system, 0j, with i = 1, . . . ,n. 
The average of the replicated partition function over the Gaussian disorder 
will be denoted by overlines. 



F = logZ= hm - (Z" - 1) . 
Now we can define an effective action by means of 



(12) 



eff 



d[(j)i] exp (-Scff [</)»]) 



with 



1=1 i=l 



i=l 



V \ - 



i=l 



(13) 



(14) 



where u = —3A^. This gives us a starting point for the analytical calculation. 
The n — i> limit should be taken at the end. 

For V = the action is 0(n)-invariant. When u = the action describes 
n decoupled Ising models. We remark that u is negative and proportional to 
the dilution. It is possible to show that for a non-Gaussian distribution, terms 
associated with higher connected momenta of the distribution appear in the 
effective action. The s-momentum couples to (p"^^ thus, if s > 2, is irrelevant 
in four dimensions and can be neglected. In our numerical simulation a site 
is occupied with probability p, so = p{l — p). 

The action ([T^) was studied in ref. [0 by using PRG techniques. Consid- 
ering a differential dilatation, the following equations are obtained: 



dr 
dlog6 



2r + 4:Kdi2u + 3v){l 
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dv 



ev - UKdvi^u + 3v) , 



(15) 



dlog6 
du 



eu - 8Kdu{4u + 3v) 



dlogfe 



where e = A — d, d being the dimension, Kd is a constant that depends on 
the dimension, and b is the Renormahzation Group (RG) scahng factor. In 
the previous formulas third order terms in u, v have been neglected, and the 
n — > limit has been considered. In the following we shall set e = 0, as we 
study a four dimensional problem. 

The initial conditions of the system, typically verify \uo\ -C Vq. Therefore, 
the RG evolution, driven by eqs. (plSj), will present two interesting regimes 

1. A transient regime in which —u ~ v"^^'^. We find v{b) ~ l/logfe or 
equivalently v{L) ~ 1/logL. 

2. An asymptotic regime reached by following the RG evolution until the 
equation 4u + 3v = 0{u^) finally holds. We obtain, including the next 
term in the perturbative expansion, v?{b),v'^{b) oc l/log6, or 1/logL. 

Thus, we expect a crossover between the pure situation where the relevant 
coupling, goes to zero as v{L) ~ 1/logL and the disordered one where 
u and V are similar in magnitude and v{L) ~ l/i/logL, i.e. the relevant 
coupling {u or v) goes to zero more slowly than in the initial regime. 

Defining t = r — AK^u, eqs. (|T5|) in the asymptotic regime reduce to 



dt 



2t + SK^ut , 



dlog6 
du 



1696 




(16) 



dlogfo 
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For large 6, the solutions are 




(17) 
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where to is an integration constant and in the large b regime u{b) does not 
depend on the initial condition uq. 

Using these formulas it is possible to obtain the expressions for the cor- 
relation length, susceptibility and specific heat as functions of the reduced 
temperature. 

We find just one slight difference with ref. 0. The equation for the wave 
function renormalization, C(^)) is 

= -7^(m, v) = -8Kl{2u^ + Quv + 3v^) . (18) 

In the pure model ( is constant, but in the asymptotic region one obtains 
C{b) oc (log 6)^/^^^ which affects to the susceptibility. This correction is neg- 
ligible from a practical point of view, but it could be important in related 
models (for instance in spin glasses and percolation P]). So the formula for 
the susceptibility reads: 



X-to^ exp 
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log to 



1/2 



iogtor/^°' , (19) 



53' 

where we have used that X = H- In ref. the term ('^ is absent. 



3.1 Calculation of Binder Cumulants 



In this section we will calculate the values of g2 and (74, with the techniques 



introduced in ref. ||T0| for the study of finite geometries using Field Theoret- 
ical methods. 

The main idea is to expand the field 4>{x) in Fourier modes. In a finite 
geometry the biggest contribution comes from the zero mode. It can be shown 
that it has to be treated non perturbatively while this is not necessary for 
the rest of the modes . 



In our case the effective action for the zero mode, that we will denote as 



ipi is, in a L'^ volume and just at the MF critical point (i.e. 



0), 



and the partition function is 

2^esin) = I I JJd?/^i j exp (- 5*65 [^i]) . 



(20) 



(21) 



vi=l 
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In the asymptotic regime, the relation 4m + 3^; ^ is satisfied with good 
precision, and so: 




Zefiin) = — = / dV'i dA exp 



(22) 



where we have introduced a Gaussian A-integration in order to decouple the 
term {Y^^=i i^fY- K is possible to see that dimensionless ratios, like and 
do not depend on the specific value of f , thereby we have also fixed f = 4!/L'^ 
in the previous formula and in the rest of the section. 

We remark that the ratio between u and f in d < 4 is universal, because 
there, we have a (limiting) fixed ratio u/v. For instance, in four dimensions u 
and V go to zero with a limiting ratio, u/v ^ —4/3, whereas in d < 4, m and 
V go to the non trivial fixed points u* and f*, respectively, and u/v ^ u*/v*. 
In (i > 4 it is impossible to fix this ratio: it depends on the parameters 
in the Hamiltonian. It is possible to show that using the e-expansion one 
can obtain for a Binder cumulant the MF value (calculated in > 4) plus 
corrections that are proportional to real and positive powers of e [|10|- So we 
can fix the ratio between u and v to the four dimensions value and then do 
the computation directly in > 4 (i.e. avoiding the loop effects). 

We can perform the integrals on the variables 

Zeff(n) = ^ / dA e~''/'lo{Xr , (23) 
v37r J 

where 

/^(A) = J d^exp [A^2 - . (24) 

For the calculation of the cumulants we need to evaluate with the action 



([2^ ) the following averages 



{M^) ^ i^a), (25) 



For instance 

= (v^^cff)-' / dA e~''/-'hM)Io{Xr-' . (26) 
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The moments in the n — hmit are 



mi) = -L / dA 



/2(A) 

/o(A)J 



e-AV3 



(27) 



where a ^ h. 

Evaluating numerically the previous integrals we obtain 

^disordered _ 0.32455... , (28) 
^disordered _ 0.31024... . (29) 

We recall that the MF values for the moments in the pure case are [|10] 

{M'n = ^ , (30) 

and for the cumulants 

c,^'' = 0.40578... , (31) 
gr"" = 0. (32) 

Now the interpretation of the formulas ( p?] ) is clear. We have a A-model with 
action 5* = Xip'^ — ip^, with A distributed with a Gaussian weight exp(— A^/3). 
Averaging with this probability distribution the moments /2m(A)//o(A) of 
these A-models we obtain the right Binder cumulants for the diluted one. 
Obviously when the mass term is zero we recover the MF result (|30|) for the 
pure model. 

3.2 FSS in the diluted model 

The scaling of the singular part of the free energy in the presence of a mag- 
netic field, /lo, is 

/sing ( ro, Wo, ^^0, /io, Y ) = &~Vsing \ r{h) ,u{h) ,v{h) , h{b), y] , (33) 



where we have introduced a new coupling, the system size L, which scales 
trivially with a RG transformation. As usually the magnetic field verifies : 

dlog h{b) -._ M'^,v) .^.N 

dlog6 2 2 ■ ^ ' 
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In the asymptotic regime, the solution of (0) is 



h{b) = /io&^(log6)2i2 



(35) 



Performing a RG transformation with b = L, we keep just one degree of 
freedom (see ref. for more details). The free energy of this system in the 
asymptotic regime is 



f{r',u',h',L = l)^\og j X 

I n n / / " 



i=l 



.1=1 



3! 

i=l 



(36) 



We re-scale the 0j variables by means of (f)[ = u^^^cpi. The free energy can be 
written as 

f{r\u\h\L = l) = f'' ^ ^' 



(37) 



obtaining finally 



/sing ro, uo, Wo, ho, — ] = L ^ f 



r{L) 



h{L) 



(38) 



We remark that the u variable is a dangerous (marginally) irrelevant 
variable fl^, [T^, and we need to do with care all the analytical steps (it is 
not correct to substitute u for its asymptotically value, -u = 0, because the 
free energy depends on inverse powers of u) . 

As / is an analytical function, the n — limit can be taken in the RHS of 
([38| ) substituting r, /i, u by their limiting values. For clarity, in the following 
we will omit the L dependence in the functions r, u and h. 

To compute the thermodynamical quantities in the critical region one just 
need to take the appropriate derivatives of /sing- It will prove convenient to 
keep in mind eqs. ([1^) and (|33D, and the following Taylor expansion (which 
depends on the relations r = t + 4K^u and the fact that t{L) = whenever 
to = 0) 



to=0 



d^f (0,1,0) + 0{u'/') 



(39) 
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where di is the partial derivative with respect to the i-ih. argument. 

As we are interested in the behavior with the lattice size just at the infinite 
volume critical temperature, the susceptibility can be written as 



smg 



ho=to=0 



dh 
dhn 



2 Q2 



ho=to=0 



L\\ogL)-^ 



1 I 

106 



(40) 



The specific heat can be computed analogously 



C oc 



dtl 



ho=to=0 



ho=to=0 



~ (logL)^e-Vili°s^l . 

At zero magnetic field, the correlation length scales as 

^(ro, Mo, 1/ L) l^^^o = L^{r, u, 1) ^^^q , 

where ^{r,u, 1) must be evaluated with the free energy 
the mass squared term is 



(41) 



and so 



(e(r,n,l)L o) 



^(ro,Uo, l/L) oc 



u 



1/2 



OC u 



1/2 



to=0 



L 



u 



1/4 



L(logL)> 



(42) 

Consequently, 
(43) 

(44) 



Finally we can also compute the shift of the apparent critical temperature. 
It can be defined as the temperature where the susceptibility (or specific heat) 
measured in a finite volume shows a maximum. Using the formula (^0|) for 
the susceptibility without imposing the constraint to = we obtain 



XocL^(logL)4+w93V(r/Mi/',l,0) . 



(45) 



The maximum of x as a function of L and t is not just at to = 0, but it 
is fixed by the condition 



(46) 



i.e. the function 1, 0) has a maximum at x 
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As if: oc Tc(oo) — Tc{L), it follows that 



Te(oo) - Te(L) oc {\og 1)-^^^^ . 



(47) 



To finish this section we will report, for completeness, the finite size for- 
mulas for the same observables in the pure case |1l4| : 



^ oc L(logL)3 , 

X oc L2(logL)2 , 

C oc (logL)3 , 

T,(oo)-Te(L) oc L~\\ogL)~i 



(48) 



The latter expressions can also be obtained with the method described 
above. 



4 Numerical Methods 

The choice of the Monte Carlo (MC) update algorithm should be carefully 
considered. The Wolff single cluster method |]T5[ is the best choice for the 
pure model. However in the diluted case, small isolated clusters of spins are 
very likely to appear. Those clusters are scarcely visited by the Wolff method. 
Therefore, we have complemented the updating method with a Metropolis 
sweep per measure. We have checked that this algorithm thermalizes appro- 
priately the configurations for p > 0.5, by comparing the numerical results 
from cold and hot starts. However, for p < 0.4, equilibration by this method 
becomes really hard to achieve. This is due to the presence of intermediate 
size clusters almost isolated from the percolating one. We have then turned 
to the Swendsen-Wang (SW) algorithm which guarantees that all-sized 
spin-clusters are considered. In this way, we find complete agreement be- 
tween hot and cold starts. We have also compared the results of a SW and a 
single-cluster simulation at p = 0.5 on our largest lattice, finding compatible 
results. However, to get a statistically-independent new configuration takes 
significantly longer (in CPU time) with the SW algorithm. 

For every observable, one first averages in the sample, then averages be- 
tween different samples. Therefore, a crucial point is how long each sample 
simulation should be. Assuming full statistical independence between dif- 
ferent measures (quite possible with a cluster method), and also between 
measures taken in different samples, the variance of such a mean is 
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where Ng is the number of samples generated and Nj is the number of mea- 
sures in each Ising-model simulation. The variance between samples of the 
thermal average of our observable is cr| and, finally, a] is the average of the 
variances in each sample. On the other hand, the computational effort is 
roughly proportional to NsNj, as the computer usually spends a fixed frac- 
tion of the time measuring. It is then clear that the optimum value of Nj 
cannot be much bigger that crf/cTg. One could even be tempted to mea- 
sure just once by sample. However, we shall come back to this point when 
discussing re-weighting methods. 



4.1 Derivatives and Re- weighting Methods 

The critical curve slope (see fig. |I]) changes quite abruptly. It is nearly vertical 
for large p and almost horizontal close to the percolation threshold. It is 
therefore wise to choose a re-weighting method such that we may extrapolate 
to different P values close to the pure model, but to different p values at very 
strong dilution. Let us then comment both extrapolation methods separately. 



4.1.1 /3-extrapoIation 

The energy measures allow the calculation of /3-derivatives of observables, 
and the use of the standard re-weighting methods, before the sample-average 
is performed. However, an important point should be made now, so let us 
recall how are they calculated: 



dp{0) = d^{0) = {0£-{0){£)) , (50) 

{0){(3 + A/?) = (Oe^'3^)/(e^^^) . (51) 

It is clear that both expressions are biased. For instance, the expectation 
value of eq. (|50|), when the averages are calculated with Nj measures, is 
really 



1 - 



2r 



(52) 
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Figure 1: Phase diagram of the four dimensional site-diluted Ising model. 
The points correspond to the simulated values, and the arrow indicates the 
percolation limit. 



r being the integrated autocorrelation time Jl^, which depends on the sam- 
ple. The bias for eq. ( pT] ) is also of order 2t/Nj, but terms of higher order 
in l/Nj are to be expected. 

These biases are immaterial for usual MC investigations, as the statis- 
tical error decreases with the square root of the number of measurements. 
However, in our case the bias is of order 1/Nj, while the statistical error is 
of order 1 / a/A^, which are similar in our simulations! 

The cure for this is to introduce unbiased estimators. A possible method 
would be to measure in completely independent samples, ensuring 2r = 1 for 
every sample. In this case the unbiased estimator is constructed multiplying 
by 1/(1 — l/Nj) the /^-derivatives and by other more complex functions the 
extrapolated observables. However, this would be too expensive from the 
computational point of view. The solution we find is to work with r > 1, 
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repeat the calculations with different values of Nj and extrapolate Nj ^ oo. 

Specifically, in each sample, we calculate the derivative with the full MC 
history, obtaining a number yi, the bias being proportional to 2t/Nj. We 
then consider two contiguous halves of the MC history, repeat the calculation 
for each one and average the final result. This value, 1/2, has a bias that goes 
as At/Nj. The next term, y^, is obtained with four quarters with a bias 
proportional to 8t/Ni. 

The linear extrapolation is 

Vl = 2yi - y2 , (53) 

and the quadratic one 

8 1 
Vq = - 2l/2 + • (54) 

We then average yL and yg for all the samples, checking that the difference 
is negligible compared with the statistical error. In fact, we have found that 
the slope in the {y, 1/Nj) plane is a very clean and easy measure of r, which 
we have used to achieve statistical independence between different measures. 
See fig. 1^ for an example. 

We proceed analogously with the extrapolations of observables or their 
derivatives. 

Another approach to eliminate the biases is to split the measures in sta- 
tistically independent sets and multiply the average of some operator in the 
first set by the average of another operator in the second set. In some cases 
(for instance derivatives extrapolated at different couplings) it could be nec- 
essary to work with more than two independent sets. As in practice one 
have contiguous MC measures, the statistical independence of the measures 
becomes involved when splitting in several sets. 

The comparison between the quadratic extrapolation and the linear one 
makes it possible to monitor this short MC history effects. In some cases 
due to the necessity of a large extrapolation in the coupling, we have found 
differences between y^ and yg around one half of the statistical error; then 
we have repeated the simulation at a nearer point. 



4.1.2 ^-extrapolation 

In addition to the standard /^-extrapolation , it is also possible to extrap- 
olate the mean values obtained at a given dilution probability, p, to a close 



15 



k 




Figure 2: In the upper part we show the sample-averaged yi value for 
taking one out of each k measures, in a L = 48 lattice at p = 0.5, with 
r ^ 8. As both r and Ni get divided by /c, we obtain an stable value until 
k = 10. In the lower part, we plot the sample-averaged i/i, y2 and values, 
as a function of the inverse number of MC measures, in a L = 32 lattice, at 
p = 0.5, with r ~ 0.8. The linear behavior is apparent. 



one p' (see [0). Let us simply recall that the probability of finding a occu- 
pation number q, when filling sites with probability p, is binomial. Therefore 
the re-weighting is calculated from a set of Ns mean values of an observable 
O and the actual density of the configuration gj)}: 

(0>(f'.« = ]v;I^(^) [t^) <0>.(«- (55) 

The visible region is of course constrained to the variance of the binomial 
distribution p(l — p)/V. Fortunately, it has been enough for us. 

Using equation (^) p-derivatives of observables can also be computed. 
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but statistical errors are about eight times bigger than for /^-derivatives. 
Therefore, our choice is to study p-extrapolated /3-derivatives, where all 
the above comments for the bias in the derivative are in order, but the re- 
weighting is unbiased. 



4.2 Numerical FSS techniques 

As already stated in the introduction, we consider the possibility of finding 
a non-Gaussian fixed point in four dimensions, where hyperscaling relations 
were fulfilled. In such the usual FSS ansatz is expected to hold. 

However, the PRG analysis rather suggests the presence of logarithmic cor- 
rections to the Gaussian behavior. Therefore, we should also consider the 
modifications in the FSS ansatz induced by the logarithmic corrections. In 
this section, we shall first remind how critical exponents are measured (see 



refs. 0, [19[ for similar calculations), then we shall show how to deal with 
logarithmic corrections. 

When hyperscaling holds, a very accurate way of measuring critical ex- 
ponents involves a form of the FSS ansatz where everything is directly mea- 
surable on a lattice: 

Om(3,p)=L^o/'^ {FoiaL,P,p)/L) + OiL--)) , (56) 

where a critical behavior t~^° is expected for the operator O, u is the uni- 
versal scaling-corrections exponent, and Fq is a (smooth) scaling function. 
Notice that terms of order ^^^^ are dropped from eq. (|56|), so we are deep 
within the scaling region. From a Renormalization Group point of view, a; 
corresponds to the leading irrelevant operator. 

In order to calculate the critical exponents, we study the quotient of 
0{sL) and 0{L), defined as 

Qo = 0{sL,P,p)/0{L,(3,p) . (57) 

Measuring at a value of the couplings where the quotient for the correlation- 
length is s the scaling function may be eliminated and we obtain: 

Qo\Q^=, = s^o/^ + 0{L-n . (58) 

If there are logarithmic corrections to hyperscaling, the critical behavior 
of the operator O is modified to 

0(L,/5„pJcx W^(logL)'5o . (59) 
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We remark that the scahng variable is again ^{L, (3,p)/ L. This is clear from 
formulas (|43|) and (pSl): the scaling variable is rjv}!'^ that is equal to 

{E,/LY, and so one readily obtains that the scaling variable is ^/L without 
logarithmic corrections. 

For the susceptibility exponent we find from eq. (0) after some algebra 



logg 



X\Q^=s 



log s 



+ 



\ogL 



(60) 



Thus, the corrections to the rj exponent are proportional to 1/ log L as in the 
pure case. 

To estimate the logarithmic correction to the critical behavior of the (3- 
derivative of the correlation length we start from eqs. (|42D and (^3|) and 

get 

i = i{r,,uoML)=L-^ . (61) 



Taking the to derivative and using eqs. ( PTj) we obtain 

3 



dfs^^L' (logL) 



1/4 



- J exp 



3 logL 
53 



(62) 



It is easy to check that 



V log s 



+ 



VlogX 



(63) 



where the correction 0(l/v^logL) arises from the exponential term in eq. (|6^). 

For the lattice sizes simulated, we expect some dependence on the dilution 
in the coefficient of the 1 / ylogX term. In the initial regime the corrections 
are proportional to 1/logL (the pure model). Until the system forget the 
initial conditions (in particular the dependence on the dilution) the coefficient 
of the 1 / i/log L term could change. 



5 Numerical Results 

The lattice sizes that we have studied have been L = 8, 12, 16, 24 and 32. 
We have generated iV5=10,000 samples, for each lattice size, at dilution val- 
ues p = 0.8, 0.65, 0.5, 0.4, 0.3. In each sample, we measure Ni = 100 times 
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after equilibration. The number of clusters traced (or SW updates) between 
measures have been chosen to yield 2r ^ 1 (see eq. ([52|)). The pure {p = 1) 
model has also been studied for L = 8, 12, 16, 24, 32, 48 and 64, as a contrast 
of the disorder-induced effects. 

We shall present our numerical results in two steps. First we shall consider 
the conventional FSS analysis (i.e. assuming hj^erscaling), finding that the 
percolation scenario is extremely unlikely. 

After that, we shall look for hyperscaling violations in the data. We shall 
find that they can be measured, and are indeed of the same order as predicted 
by PRC. 



5.1 Assuming Hyperscaling 

To measure the critical exponents, we use the so-called quotients method, 
which allows for a great statistical accuracy 0, The starting point is eq. 
([5S|). We have first approximately located the point where 

then we have used re- weighting techniques to fine-tune the condition (|5^. 
For p = 1.0, 0.8, 0.65 and 0.5, we have used, /^-extrapolation. Therefore, the 
critical p is fixed, (3 being the tunable parameter. For lower dilution values, 
we have rather used p-extrapolation, so we have first approximately located 



the f3 value, for which condition (64) holds at p ~ 0.3, 0.4. Next we fix /3, 
fine-tuning p afterwards, so the critical values differ from p = 0.3, 0.4 by an 
amount of less than 1%. However, in tables and figures, we shall refer to 
them as p = 0.4 and 0.3 for brevity. 

Eg . (|58|) applied to the operators djs^ and yields respectively the expo- 
nents 1 + 1/z/ and 2 — 1]. The numerical results are shown in tables |I| and 
|[ For the u exponent we find, instead of a stable value, a monotonically 
decreasing one. For rj, such an evolution with growing L is found, but it is 
clearly weaker. Therefore, an infinite volume extrapolation is called for. If 
hyperscaling holds, we expect finite-volume scaling-corrections proportional 



to L~^. As 00 = 1.13(10) in the percolation [|T9|, in the last row of both 
tables and H we include an infinite- volume extrapolation with u; = 1. This 
fit is shown in figure ^. 

It is clear that the percolation scenario, u = 0.686(2) and t] = —0.094(3) 



19|, can be ruled out. Moreover, the possibility of a different fixed point 
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L 


p= 1.0 




p = 0.65 


j5==0.5 


p~ 0.4 


p~0.3 


8 


0.5119(9) 


0.5175(11) 


0.5308(13) 


0.5482(16) 


0.5604(15) 


0.5700(26) 


12 


0.5074(18) 


0.5154(11) 


0.5270(13) 


0.5428(19) 


0.5532(19) 


0.5647(22) 


16 


0.5066(10) 


0.5142(13) 


0.5251(12) 


0.5412(19) 


0.5478(18) 


0.5583(26) 


24 


0.5067(18) 












32 


0.5039(17) 












oo 


0.5019(14) 


0.5110(25) 


0.5194(26) 


0.534(4) 


0.536(4) 


0.549(5) 



Table 1: The u exponent for {L,2L) pairs at the different dilutions. 



L 


p^ 1.0 


j3 = 0.8 


p = 0.65 


p^O.5 


p~ 0.4 


p~0.3 


8 


-0.02644(14) 


-0.0161(12) 


-0.0124(11) 


-0.0055(19) 


0.002(8) 


-0.003(4) 


12 


-0.02132(26) 


-0.0138(12) 


-0.0089(13) 


-0.0051(17) 


0.000(9) 


-0.0043(21) 


16 


-0.01656(12) 


-0.0130(12) 


-0.0053(10) 


-0.0073(17) 


0.000(11) 


-0.006(7) 


24 


-0.01259(29) 












32 


-0.01445(18) 












oo 


-0.0085(18) 


-0.0097(25) 


0.0013(22) 


-0.008(4) 


-0.002(20) 


-0.007(9) 



Table 2: The rj exponent for (L, 2L) pairs at different dilutions. 



neither Gaussian nor with the percolation critical exponents requires a fairly 
exotic FSS behavior. We do not find this possibility likely. 

At this point, one could claim that this model presents weak universality, 
as recently proposed in the two dimensional version of the model That 
is, the exponent 1] is constant over the critical line, while u is continuously 
varying. However, a much less spectacular, but more likely interpretation 
will be given in the next subsection. 



5.2 The quest for logarithms 

As PRG predicts logarithmic corrections to the MF behavior, we should 
expect scaling-corrections of order 1 / log L (for the 1] and u exponents in the 
pure model and only for the rj exponent in the diluted case) and 1 / -y/log L 
(for the u exponent in the diluted one): both are of the same order for the 
lattices that we can afford!. In figures ^ and ^ we show that the deviation 
from MF can indeed be accounted for by logarithmic corrections. We have 
found similar results in two dimensions ||21| . While finishing these papers, the 
same conclusion has been independently drawn in a transfer-matrix study of 
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0.50 




0.00 



0.05 



0.10 



0.15 



l/L 



Figure 3: Exponent v for different lattice pairs of type (L, 2L). The lines 
correspond to linear fits. 



the bond-diluted two dimensional model p2 . 

As we have seen, hyperscaling does not seem to hold. Indeed, the PRG 
predict logarithmic violations. In this section, we try to identify them by 
starting from a naive point of view. That is, we shall first locate the critical 
point as if the usual FSSA were correct, and then we shall study there the 
scaling of physical quantities, looking for deviations from the pure-power law. 

A very accurate way of computing (3c (or pc) is to fix a value of g^, 
measuring in what [3 the function g^^L,^) is equal to the fixed previous 
value p3. In absence of logarithmic corrections one expects that the shift in 



P behaves as In figure ^ we plot these quantities for three values of g^ 

at each p, while the results of a quadratic fit, for the infinite volume critical 
couplings are found in table 0. We point out that a logarithmic correction 
as the one computed in ref. |]I4| (as eq. (|49|) shows) for the pure model or 
eq. ( ^71) for the diluted one, does not change the fit results. 
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0.0 



0.2 



0.4 0.6 

l/(logL)'/^ 



0.8 



Figure 4: The v exponent obtained from (L, 2L) pairs. The sohd hnes are 
hnear fits constrained to be z/ = 0.5 in the L oo hmit. 

In order to minimize the systematic errors, those values have been ob- 
tained choosing a such that the hnear coefficient in vanishes. We 
have also computed the extrapolations for a wide range of values to check 
the amplitude of the change. We observe that, within one standard deviation 
in the extrapolated value, (74 can be changed in the interval [0.3, 0.8] in the 
best case [p = 0.8) and in the interval [0.4,0.6] in the worst one {p = 0.5). 
Although we think that systematic errors from this source are negligible, a 
more conservative attitude could be to duplicate the statistical error bars in 
table 1^. 

Now, we can check hyperscaling violations of the form: 

e(i:,/?,)cxL(logL)'« . (65) 

In figure 0, we plot log(^/L) as a function of log(logL), the slope being 
directly 6^. We can see that a good linear behavior is obtained. The fitted 6^- 
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0.02 



0.00 



-0.02 




0.02 —73 = 



0.00 



-p = 0.65 
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, , , , F 



1/IopL 



Figure 5: The 77 exponent for different lattice sizes. The fit is enforced to 
yield 77 = for L — > 00. 



values are reasonably close to the theoretical prediction, but a growing trend 
with the dilution is self-evident. A naive (wrong) explanation is that this is 
an effect coming from the vicinity of the percolation critical point. Indeed, we 
can repeat the previous fit in the pure percolation for similarly-sized lattices. 
Fitting from L = 8 we obtain = 0.0927(9) with x^/d.o.f. = 97.9 and from 
L = 16, the fit parameters are = 0.066(2) with x^/d.o.f. = 6.90. So a 
linear behavior is ruled out for this model. However, this is not surprising, 
as hyperscaling is known to hold for percolation in four dimensions. 

Another interesting observable is the specific heat, which for the pure 
model is expected to diverge as C oc (logL)^*^^, while in the diluted case it 
is expected to remain bounded. Both predictions can be tested. The fitted 
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Figure 6: Computation of (3^ or looking at fixed values of (74 for the different 
values of the dilution. The lines are quadratic fits. 

values for the logarithmic divergence exponent, 6c ■, are 



V = 


1.0 


Sc 


= 0.399(4 + 22) 


p = 


0.8 


Sc 


= 0.304(7+ 13) 


p = 


0.65 


Sc 


= 0.184(8 + 15) 


p = 


0.5 


Sc 


= 0.095(6 + 9) 


p = 


0.4 


Sc 


= 0.084(5 + 6) 


p = 


0.3 


Sc 


= 0.073(8 + 7) 



(66) 



On the values of the Binder cumulants at the infinite volume critical 
temperature we can see directly the (logarithmic) corrections to the scaling Q. 
In the pure case we have obtained ~ 0.51 for the largest lattice at the 
critical point: we remark that the asymptotic value for the pure model is 

^ Do not confuse with the multipUcative logarithmic corrections to the critical law that 
we have studied above. 
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Pc 


Pc 


1.0 


0.149695(1) 


0.8 


0.188864(3) 


0.65 


0.235049(8) 


0.50 


0.317368(19) 


0.398806(18) 


0.42 


0.30110(4) 


0.633 



Table 3: Critical couplings as obtained from a quadratic fit of Pc{L, g^). The 
error bars are purely statistical. 



5^4 ~ 0.406. It is possible to show using that the leading correction to 
scaling term for goes as 1/ log L in both the pure model and the diluted 
one. Our pure (?4-data are compatible with a limiting value of 0.406 modified 
by 1/logL corrections. 

In the diluted model, we have obtained data for the (74 cumulant around 
0.5 (see fig. ^) which is quite larger than the predicted value g^ ~ 0.32 (see 
eq. (|29|) ). Again, this numerical discrepancy can be understood taking into 
account corrections like 1/logL. The same comments hold for g2. 

In figure ^jwe have shown g2 against g^^ for some values of the dilution. We 
have also plotted the theoretical value for the pair {g2,gi)- This figure must 
be interpreted having in mind the previously cited logarithmic corrections to 
the scaling. It is clear that the complete numerical characterization of these 
cumulants needs further numerical work. 

Finally we have studied the probability distribution of the observable 
Ai^ (i.e. we collect the histogram of the values of Ai^ at every measure 
independently of the sample). We have found that for large values of Ai^ 
the data follow a law P(A^^) oc exp(— c(L)(A^^)^). Using the results of 
section (|3.1|) it is possible to show that the theoretical prediction for the 
probability distribution is a Gaussian with c{L) oc u{L)L^. The coefficient, 
c(L), computed numerically follows very well a law L^/^/[ogL in perfect 
agreement with the theoretical prediction. 
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-0.6 



p = 0.5 

6f=0.189(7+12) 



p = 0.4 

5f=0.241(8+12) 




5,=0. 181(6+14) 
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0.8 1.0 1.2 



0.8 1.0 1.2 



log(logL) 



Figure 7: log(^/L) in the infinite volume critical point, as a function of 
log(logL), for several dilutions. The fitted values are also displayed. The 
first term in the error is statistical, while the second corresponds to the error 
in the critical coupling. The three lines correspond to data at the critical 
point and one standard deviation apart at each side. 



6 Conclusions 

We have performed a Monte Carlo simulation of the site-diluted Ising model 
in four dimensions, for several values of the dilution, and in a wide range of 
lattice sizes. The use of a finite-size scaling analysis allows us to consider 
big lattices just at the critical point. To gain accuracy we have repeated the 
simulations for many different hole configurations. 

As a first stage, we have measured with great precision the critical expo- 
nents under the hypothesis of hyperscaling. The value we obtain for the i' 
exponent changes along the critical line, ruling out the possibility of a single 
non-Gaussian fixed point. 
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Figure 8: g/^ versus §2 for three values of the dilution. The symbol sizes 
increase with the lattice size (L = 8, 12, 16, 24, 32). The PRG prediction is 
also plotted. 

Using Pcrturbative Renormalization Group techniques, we have com- 
puted the scaling formulas for the diluted model, obtaining specific loga- 
rithmic corrections to the Gaussian behavior. 

We have re-analyzed our numerical data finding that they agree with 
a Gaussian scenario with logarithmic corrections to hyperscaling. The pure 
model has also been considered and its corresponding scaling formulas, checked. 

Although the nature of the logarithmic corrections hardly allows to per- 
form precise fits to the predicted functional forms, we have found a reasonable 
agreement. 
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